from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import geopandas as gpd
import altair as alt
%matplotlib inline
pd.options.display.max_columns = 999
plt.rcParams['figure.figsize'] = (10,6)
Oct 29, 2020

Three parts:
import osmnx as ox
Make sure the version >= 0.16 for this notebook to work!
print(ox.__version__)
Haven't we already done this?
We've used hvplot, holoviews, datashader, etc. to create interactive map-based visualizations
Pros
Cons
OSMnx leverages Folium under the hood to make interactive graphs of street networks!
Key function: ox.plot_graph_folium will make an interactive map of the graph object
Load the street network around City Hall
G = ox.graph_from_address('City Hall, Philadelphia, USA',
dist=1500,
network_type='drive')
# plot the street network with folium
graph_map = ox.plot_graph_folium(G,
popup_attribute='name',
edge_width=2)
ox.plot_graph_folium?
And now save the map object and load it into the Jupyter notebook
from IPython.display import IFrame # loads HTML files
filepath = 'graph.html'
graph_map.save(filepath)
IFrame(filepath, width=600, height=500)
Folium map objects are supposed to render automatically in the Jupyter notebook, so if you output a Folium map from a notebook cell, it will render.
However, there a lot of times when the map won't show up properly, especially if the map has a large amount of data. Saving the file locally and loading it to the notebook via an IFrame will always work.
type(graph_map)
graph_map
Let's calculate the shortest route between the Art Museum and the Liberty Bell.
See last lecture (Lecture 9A) for a guide!
Use OSMnx to download all amenities in Philadelphia of type "tourism".
ox.geometries_from_place() can download OSM features with a specific tag
You should notice we have the building footprint for the Art Museum (a polygon geometry) and the point location for the Liberty Bell.
.squeeze() function can be useful for converting to a Series object from a DataFrame of length 1
For the Art Museum geometry, we can use the .geometry.centroid attribute to calculate the centroid of the building footprint.
Remember: we need the coordinates in the order of (latitude, longitude)
Use the street network graph around City Hall and the ox.get_nearest_node() function to find the starting and ending nodes for the trip.
networkx to find the shortest path between nodes¶The nx.shortest_path() will do the calculation for you!
import networkx as nx
Now, we can overlay the shortest route between two nodes on the folium map.
Key function: use ox.plot_route_folium to plot the route.
# plot the route with folium
route_map = ox.plot_route_folium(G, route)
filepath = 'route.html'
route_map.save(filepath)
IFrame(filepath, width=600, height=500)
We can also add the underlying street network graph
# plot the route with folium on top of the previously created graph_map
route_graph_map = ox.plot_route_folium(G, route, route_map=graph_map)
# save as html file then display map as an iframe
filepath = 'route_graph.html'
route_graph_map.save(filepath)
IFrame(filepath, width=600, height=500)
Note the Leaflet annotation in the bottom right corner of the maps...
Things we'll cover:
import folium
Key function: folium.Map
Some key ones:
Let's take a look at the help message:
folium.Map?
# let's center the map on Philadelphia
m = folium.Map(
location=[39.99, -75.13],
zoom_start=11
)
m
m = folium.Map(
location=[39.99, -75.13],
zoom_start=11,
tiles='stamenwatercolor'
)
m
Let's try out the ESRI World Map:
Important: for custom tile providers, you need to specify the attribution too!
tile_url = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer/tile/{z}/{y}/{x}'
attr = 'Tiles © Esri — Source: Esri, DeLorme, NAVTEQ, USGS, Intermap, iPC, NRCAN, Esri Japan, METI, Esri China (Hong Kong), Esri (Thailand), TomTom, 2012'
m = folium.Map(
location=[39.99, -75.13],
zoom_start=11,
tiles=tile_url,
attr=attr
)
m
Key function: folium.GeoJson
Key parameters:
style_function: set the default style of the featureshighlight_function: set the style when the mouse hovers over the features tooltip: add a tooltip when hovering over a featurefolium.GeoJson?
folium.GeoJsonTooltip?
# Load neighborhoods from GitHub
url = "https://github.com/azavea/geo-data/raw/master/Neighborhoods_Philadelphia/Neighborhoods_Philadelphia.geojson"
hoods = gpd.read_file(url).rename(columns={"mapname": "neighborhood"})
# Load ZIP codes from Open Data Philly
zip_url = "http://data.phl.opendata.arcgis.com/datasets/b54ec5210cee41c3a884c9086f7af1be_0.geojson"
zip_codes = gpd.read_file(zip_url).rename(columns={"CODE":"ZIP Code"})
ax = ox.project_gdf(hoods).plot(fc="lightblue", ec="gray")
ax.set_axis_off()
ax = ox.project_gdf(zip_codes).plot(fc="lightblue", ec="gray")
ax.set_axis_off()
Define functions to set the styles:
def get_zip_code_style(feature):
"""Return a style dict."""
return {"weight": 2, "color": "white"}
def get_neighborhood_style(feature):
"""Return a style dict."""
return {"weight": 2, "color": "lightblue", "fillOpacity": 0.1}
def get_highlighted_style(feature):
"""Return a style dict when highlighting a feature."""
return {"weight": 2, "color": "red"}
# Create the map
m = folium.Map(
location=[39.99, -75.13],
tiles='Cartodb dark_matter',
zoom_start=11
)
# Add the ZIP Codes GeoJson to the map
folium.GeoJson(
zip_codes.to_crs(epsg=4326).to_json(), # IMPORTANT: make sure CRS is lat/lng (EPSG=4326)
name='Philadelphia ZIP_codes',
style_function=get_zip_code_style,
highlight_function=get_highlighted_style,
tooltip=folium.GeoJsonTooltip(['ZIP Code'])
).add_to(m)
# Add a SECOND layer for neighborhoods
folium.GeoJson(
hoods.to_crs(epsg=4326).to_json(), # IMPORTANT: make sure CRS is lat/lng (EPSG=4326)
name='Neighborhoods',
style_function=get_neighborhood_style,
tooltip=folium.GeoJsonTooltip(['neighborhood'])
).add_to(m)
# Also add option to toggle layers
folium.LayerControl().add_to(m)
m
.to_json()LayerControl to toggle different layers on the mapfolium.GeoJsonTooltipOverlay GeoJSON features on an interactive map, colored by a specific data variable
In the interest of time, I've used cenpy to download data for internet availability from the 2018 5-year ACS, but a good at-home exercise is to try to replicate my work.
census_data = pd.read_csv("./data/internet_avail_census.csv", dtype={"geoid": str})
# Remove counties with no households
valid = census_data['universe'] > 0
census_data = census_data.loc[valid]
# Calculate the percent without internet
census_data['percent_no_internet'] = census_data['no_internet'] / census_data['universe']
Load counties from the data folder as well:
counties = gpd.read_file("./data/us-counties-10m.geojson")
folium.Choropleth?
Steps:
key_on=), using the GeoJSON "features.properties." syntaxcolumns= keywordm = folium.Map(location=[40, -98], zoom_start=4)
# Convert the counties geometries into GeoJSON
counties_geojson = counties.to_crs(epsg=4326).to_json()
folium.Choropleth(
geo_data=counties_geojson, # Pass in GeoJSON data for counties
name='choropleth',
data=census_data, # the census data
columns=["geoid", 'percent_no_internet'], # First column must be the key, second the values
key_on="feature.properties.id", # Key to match on in the geometries
fill_color='RdPu', # any ColorBrewer name will work here
fill_opacity=0.7,
line_opacity=1,
line_weight=0.5,
legend_name='Households without Internet (%)'
).add_to(m)
m
census_joined = pd.merge(
counties[["id", "geometry"]], census_data, left_on="id", right_on="geoid"
)
census_joined.head()
We will use a matplotlib colorbar and it requires data to be between 0 and 1
x = census_joined['percent_no_internet']
census_joined['percent_no_internet_normalized'] = (x - x.min()) / (x.max() - x.min())
import matplotlib.colors as mcolors
# use a red-purple colorbrewer color scheme
cmap = plt.get_cmap('RdPu')
cmap(0.5)
mcolors.rgb2hex(cmap(0))
mcolors.rgb2hex(cmap(1.0))
def get_style(feature):
# get the data value from the feature
value = feature['properties']['percent_no_internet_normalized']
# evaluate the color map
# NOTE: value must between 0 and 1
rgb_color = cmap(value) # this is a RGB tuple
# convert to hex string
color = mcolors.rgb2hex(rgb_color)
# return the style dictionary
return {'weight': 0.25, 'color': color, 'fillColor': color, "fillOpacity": 0.75}
Be sure to only keep the data columns we'll actually use!
needed_cols = ['NAME', 'percent_no_internet', 'percent_no_internet_normalized', 'geometry']
census_json = census_joined[needed_cols].to_json()
# initialize the map
m = folium.Map(location=[40, -98], zoom_start=4)
# add the GeoJson to the map
folium.GeoJson(
census_json,
name='choropleth',
style_function=get_style,
highlight_function=get_highlighted_style,
tooltip=folium.GeoJsonTooltip(['NAME', 'percent_no_internet'])
).add_to(m)
folium.LayerControl().add_to(m)
# avoid a rendering bug by saving as HTML and re-loading
m.save('percent_no_internet.html')
The hard way is harder, but we have a tooltip and highlight interactivity!
IFrame('percent_no_internet.html', width=800, height=500)
Try to replicate the above interactive map exactly (minus the background tiles). This includes:
Altair's syntax is similar to the folium.Choropleth syntax — you should pass the counties GeoDataFrame to the alt.Chart() and then use the transform_lookup() to merge those geometries to the census data and pull in the census data column we need ("percent_without_internet").
Hints
import altair as alt
One of leaflet's strengths: a rich set of open-source plugins
https://leafletjs.com/plugins.html
Many of these are available in Folium!
from folium.plugins import HeatMap
HeatMap?
New construction requires a building permit, which we can pull from Open Data Philly.
carto2gpd package to get the dataYou can use the "current_date" variable (built in to the SQL database) to get the current date
For example, to get the shootings from the past 7 days, we could do:
SELECT * FROM permits WHERE permitissuedate >= current_date - 7
To select new construction permits, you can use the "permitdescription" field, there are two relevant categories:
IN command (documentation)import carto2gpd
Some shooting incidents don't have locations — use the .geometry.notnull() function to trim the data frame to those incidents with valid geometries.
Note: We need an array of (latitude, longitude) pairs. Be careful about the order!
The HeatMap takes the list of coordinates (first column lat and second column lng)